% plot the g/f curve for an izhikevich neuron.
%function [slope_intercept plane asave es is]=gfcurve

% names{1}='Pyr old';
% names{2}='Pyr 2008';
% names{3}='FS 2008';
% names{4}='TC 2008';

names{1}='Adaptive Exponential pyr.';
names{2}='Adaptive Exponential FS';

n=1;
for neuron = 1
    figure(1);
    j=1;
    es=0:.2:15;
    is=0:.2:25;    
    for e = es
        k=1;
        for i = is
%            a(j, k, neuron)=izhikevich(e,i,neuron,true);
            a(j, k, n)=adaptexp(e,i,neuron,true);
            %a(j, k, n)=adaptexp(e,i,ga_params,true);
            k=k+1;
        end 
        j=j+1;
    end

    figure(neuron+1);
    imagesc(a(:,:,n));
    title (['Conductance/firing rate parameter for neuron' names{neuron}],'FontSize',14);
    colorbar;
    xlabel('I g (nS)','FontSize',12);
    ylabel ('E g (nS)','FontSize',12);

    drawnow
    n=n+1;
end

neuron
z=a(:,:,1);

a(a==0)=10^6;[junk e]=min(a(:,:,1));
slope_intercept=polyfit(es(e),is,1);
plot(es(e),is,'.');

asave = a;

% Now determine the rate equation as a function of ge and gi by fitting
% a plane to the data.
[x y]=meshgrid(is,es);
% dilate the array of nonzero elements to include the border of zeros to
% get a better fit at the bottom end of the curve.
se2 = strel('line',3,90)
ind = find(imdilate(z>0.000001 & z<.75*max(max(z)),se2));
[a b c]=planefit(x(ind),y(ind),z(ind))
plane = [a b c];
figure
subplot(1,2,1);
surf(z)
subplot(1,2,2);
planez=a*x+b*y+c;
planez(planez<0)=0;
surf(planez);
figure;
surf(z-planez);


